Summary created by Matthias Meyer. Note there might be errors!

Contents

Some parts are also available on the following webpage.

1 Newton Polynomische Interpolation

1.1 Collocation

Collocation = All measurement points are represented by a function, for example a polynomial. The polynomial

results in a linear equation system with a degree of n+1.

1.2 Aitken-Neville Rekursionsformel

One way to solve this is by searching a polynomial formula as can be seen below. (It would also possible to generate a function by connecting the different data points with a straight line. ⇒ On a Computer this would take a lot of computational effort, since one would have a lot of if and else statements)

To get the solution of this polynomial, one has to solve the following equation system:

As one can see, this equation system gets really large and could be difficult to be solved on a microcontroller. But there exists a nice algorithm which makes solving this system easier, called Aitken-Neville recursion which divides the huge equation system in little parts.

The formula above shows how the global interpolating polynomial is combined from the partial interpolation polynomials p_{0,1,2,…,n−1}(x)and p_{1,2,…,n−1,n}(x). On this partial interpolation polynomials, one can again apply the formula until one ends up with only two datapoints.

1.3 Newton basis polynomials

Since the calculation with the Aitken-Neville recursion is quite tedious. Newton came up with another basis polynomial π_k(x) with k= 0,1,2,...,n (Aitken-Neville recursion used (1,x¹,x²,....,x^m)

(1)

Which results in the final polynomial which can be seen in Equation 2

(2)

When one now writes the equation system one sees that this system is much easier to solve:

When one also applies the Aitken-Neville recursion it even get easier and independent of the order, since one just calculates divided differences. Below one can see the calculation of the first a_k terms.

Where y is called the divided difference.

When the points have the same distance to each other, the formula gets even easier:

1.4 Example: Collocation polynomial

Given is the following dataset:

Calculate the collocation polynomial.

One can do the calculation the following way:

and therefore get the following result for the given data points:

According to Equation 2 one gets then the following result:

1.5 Further information

It does not depend on which data point one uses first and which one as last element. The resulting formula might look different (different a_k coefficients, but the last one is the same), but the result is exactly the same. Furthermore, the data points must not have the same spacing.

1.6 Problems

With a lot of data point the runge phenomenon occurs (oscillations with high frequencies and amplitudes towards the boundaries of the arguments range) To calculate the error which can be seen one can use the following formula:

Where ξ is a new data point.
can also be substituted by ⇒C=
It is a newton polynomial with one more argument. is a higher order derivative which we do not know at the moment. The formula can also be rewritten in Equation 3:

(3)

When generalizing it one can also write Equation 4

(4)

1.6.1 Error Calculation example

The model fucntion y(x) =sin(πx) has to be interpolated using the arguments x= 0;1;2 by a quadratic polynomial p(x), what is the error at the position x=

Error = y(x) −p(x) =sin(πx) −_{0+1(x−0)+1(x−0)(x−1)}

2 Chebyshev arguments

To reduce the error mentioned above to a minimum one can use a chebyshev distribution of the arguments.

2.0.1 Exercise 2 and 3 of excercise sheet 2 are important for the exam!

3 Hermite Interplation (Osculation)

When the problem of collocation, which can be found on the following post is extended by the requirement that certain given values of derivatives of order 0 up to some higher order k of the model function y must be met at some of the arguments x₀,x₁,...,x_n we end in an interpolation problem called osculation or Hermite interpolation. Furthermore note that one must use the modified newton polynomials as they can be seen in Equation 5, for the example provided there

3.1 Example: Hermite Interpolation and Error Callculation

We have railway track with the following points given:(0;0),(2;1),(4;2) and also its derivatives. Now we have to search a polynomial p₂ which goes through point one and point two and fulfils it’s derivatives. So we have the following conditions: p₂(2) =1,p₂′(2) =1,p₂′′(2) =0 and p₂′(4) =0,p₂(4) = 2,p₂′′(4) =0

This lead to the following result: p₂ =1 ·π₀ +1 ·π₁ + 0 ·π₂ −π₃ +π₄ when using the modified newton polynomials as they can be seen in Equation 5.

(5)

The problem with this method it is that it is not guarateed to find a solution for this problem, when not all derivatives are given. One then has to increase the order of the polynomial.

Error With Equation 4 the error is then:

The maximum error therefore is: ··

3.2 Example 2

Compute two fourth-degree (!) polynomials, p₁(x) and p₂(x), meeting the constraints below:

and

Moreover, the two polynomials should meet smoothly at the point (2,1) without a crinkle (with a common tangential line) To solve this problem we introduce a new variable called a which defines the first derivative at point (2,1). When a is equal in both equations the meeting of the two polynomials is very smoothly. The first scheme looks like this:

Where

since p₁ must be of order 4 se conclude that |F|= 0(!) and therefore a= 1

Now one can do the same thing for the next polynomial but a is this time known. When solving it one gets the following result:

4 Multi-variable Polynomial Interpolation

The polynomial interpolation can also be used with a Multi-variate Polynomial Interpolation. Where the polynomial is represented by Equation 6.

(6)

4.1 Example: Multi-variable polynomial Interpolation

p(x,y) p(0,0) =0,p(1,0) =1,p(0,1) = 0,p(1,1) = 0.5 Now lets calculate the first x row:

Now lets calculate the second x row:

And in step 3 we combine those two.

4.1.1 Alternative Method:

4.2 Questions

How many conditions are generally required

• for a tri-linear polynomial interpolation

  • We have three dimensions and a linear polynomial ⇒ we have 8 basis functions (2 · 2 · 2) ⇒ 8 conditions that are required:

    

• for a tri-cubic polynomial interpolation

  • We have again three dimensions and each dimension has four newton basis polynomials ⇒ 4 · 4 · 4 = 64

  • We would need to define;

    

5 Spline Interpolation

5.1 Idea

The Idea of the spline interpolation is that one does not interpolate the data with a high degree polynomial, but with multiple polynomials of lower degree. Due to that the runge phenomenon does not occur which occurs for high deggre interplations. When following this approach the transition of one spline ot the next must be considered. Normally one says that the derivative up ot an order n must be the same from one to the next spline. The drawback of this approach is that one needs a lot of storage, since one needs to store a lot of fucntions. The advantage is that it is easier to calculate, since the newton interpolation has a complexity of n² whereas a cubic spline interpolation has a complexity of n.

5.2 Cubic Spline

5.2.1 Solve problem

A spline can be described with Equation 7

(7)

5.2.2 Natrual Splining

In natural splines the energy is minimized, therefore y₀^′′= 0 =y_n^′′

• For natural spline c₀ and c_n are given by Equation 8

  (8)

• For cubic splines the a coefficients can be calculated according to Equation 9

  (9)

• The b coefficients can be calculated according to Equation 10

  (10)

• The d coefficients are given by Equation 11

  (11)

• solve equation system after c_i:

(12)

Error Calculation for cubic splines with C²

(13)

5.2.3 Formulas for the cubic clamped spline interpolation S

• b₀ =y₀^′;b_n=y_n^′
• solve equation system after c_i:

(14)

• b_i−1 = −h_i−1 (i= 1,…n− 1)
• b_n−1 = −c_n−1h_n−1 −d_n−1h_n−1² = −c_n−1h_n−1
• d_i−1 = (i= 1,…n− 1)
• Clamped splining: d_n−1 = Natural Splining: d_n−1 =−

5.2.4 Example Natural Spline

Calculate the natural cubic spline-interpolation for a sine function sin(x) in the interval [0,π] according to the points {0,,π}. Also calculate the maximum error.
From Equation 7 one knows that

Furthermore c₀ = 0 (Equation 8) since we use natural splines and a₀ =y₀ = 0,a₁ =y₁ = 1 (Equation 9). From Equation 13 one can write down the following equations:

From Equation 10 one knows that

From Equation 11 one knows that

And finally from Equation 10 that:

The error can then be estimated with Equation 13.

5.2.5 Example

The data below was generated by the sine function. In this example, the natural and clamped spline as well as the max error are calculated.

natural Spline

(15)

clamped Spline

(16)

For the two examples above (natural and clamped) give maximum estimations for the following error quantities ∣y(x) −S(x)∣, ∣y′(x) −S′(x)∣ and ∣y′′(x) −S′′(x)∣. The osculation error of a cubic spline can be calculated with ?? and the osculation error of a periodic cubic spline () with Equation 13.,

5.3 Bernstein-Bézier Splines (B-B-Splines)

The bernstein-Bézier splines should give the same result as the cubic splines mentioned in the previous chapter. The difference is that one does not get a single formula in the end, but different data points. A good explanation can be found in the following video. But first of all to understand bézier curves/spline one must be familiar with Bernstein polynomials and therefore with the binomial coefficient (see also Equation 17)

(17)

Example: (For what can the binomial coefficients be used?) What is the fourth term of (3x− 4y)⁶ (note there is no zeroth term, therefore we subtract one in the equation blew)? The result can also be read from Figure 2.

Therefore, the fourth term is:

This is much easier than to really calculate the polynomial. An explanation can also be found in the following video

5.3.1 Bernstein Polynomial

The bernstein polynomial is defined in Equation 18.

(18)

In Equation 18 one had an interval ∈ [0,1]. When our data is in another range like ∈ [a,b] one has to use Equation 19

(19)

Furthermore note that those polynomials look like one can see in Figure 3.

The polynomials relate to each other, as one can see in Equation 20.

(20)

5.3.2 Simple Bézier Curves

A simple Bézier curve is defined with Equation 21. To get the idea also have a look at Figure 4.

(21)

Where:

• i control point number
• n is the total number of control points minus one, since Points are (P₀,P₁,P₂,P_n)

5.3.3 Composite Bézier Curves

The simple Bézier curves meet at common control points. Which is a continuity condition (C⁰), but often higher (smoothness) conditions are required (C^k-smooth). This condition is only met if and only if Equation 22 is given.

(22)

Writing out Equation 22 for C¹ smoothness results in Equation 23, whereas C² smoothness results in Equation 24 where one has to know that also Equation 23 C¹ smoothness must be met.

(23)

(24)

The Bézier-curve is defined by the control points and the Bernstein polynomials. On each spline one has n+1 control points.

• spline number, the spline number there are m splines (m=(Number of given Data points-1))
• point on spline, there are n+ 1 points on the spline

5.3.4 Example: Composite Bézier Curves

The four points A= (0,0),B= (1,0),C= (2,3) and D= (2,4) are to be interpolated (joined) by composed C¹ Bernstein-Bézier splines: A and B are to be joined linearly (by a straight line), as well as C and D.
Compute the missing C¹ Bernstein-Bézier spline of minimal degree between B and C.

To solve the exercise one can use Equation 23. Where the first spline and the last one have two control points, since it is a straight line ⇒n=1. The second spline has four conditions, since tow points must be met and two derivatives, since it must be C¹ smooth. Due to that n=3 ⇒ the degree is also 3 (also called cubic).

5.3.5 Properties

• The Bézier-curves is always inside the convex hull of the data points.

• 
• If C^k smoothness for a point is required, k equations or k control points are necessary

5.3.6 Casteljau recurrence

The Casteljau recurrence is a similar idea as the neville-aitken. With this Idea a point on the Bézier curve can be calculated as a linear combination of two points on a Bézier curve of a lower degree.

Is (t) =∑ _i=0ⁿP_i,jB_i,n(t) t∈ [0,1] with the functions Q_j defined (degree: n ), the following formulas turn out:

5.3.7 Example

Let’s do the same example as in subsubsection 5.2.5. Where the following four points are given:

Therefore Q₀ = (0,0),Q₁ =,… and h_j=h=. One now has to met the following requirements:

The solution of the linear system above gives us the following points:

When we now calculate the first spline we get the same result as before.

(25)

For the second spline we get the following:

(26)

6 Linear Least-Squares approximation

6.1 Idea

Interpolation with the collocation methods often run into oscillation problems for (rather large) sets of measurement points. Furthermore, in most cases the measurements also contain some error points, which one does not want to represent in the graph. Due to that, an approximation (data points are not represented exactly any more) might be the preferred way to represent the data.

6.2 Linear Least-Squares

To find the best approximation, one must define what is a good and what is a bad approximation, which could be with the following basic functions:

Note: Error = residuals ⇒ norm of residuals
Mathematically, the minimization of the squared errors is the easiest, therefore this one is most commonly used (least square approximation)
Figure 6: Data approximated by an error curve

The approximation function can be described with a set of basis functions, which we name here g₀,g₁,…,g_m=_j=0,…,m. Note the basis functions are sometimes also called monomials. Furthermore we define the following

• N: Number of sample points
• m: degree of basic functions (m ≪ N)
• a: weighting coefficients
• ,₁,…,=_i=0,…,N: measurement points

With those variables one can create a equation in matrix notation form, as can be seen in Equation 27:

(27)

Since Equation 27 is normally overdetermined (m≪N) The error/residuals can be calculated with Equation 28 and the squared sum S of residuals with Equation 29. The goal is now to minimize S from Equation 29. This can be done with Equation 30, which is not derived in this post (for more information, search after orthogonal projection).

(28)

(29)

6.2.1 Thinking hint

Lets assume one has the following points: {1,1},{2,2} and one wants to approximate those points by a polynomial of degree zero (m= 0) ⇒g₀(x) = 1. Then one can write the term inside the square brackets as the following:

But as one can see one can not take the square of this term, therefore one has to multiply on both left sides with G^T which has no effect on the end result, since one is only interested in the minimum.

When one squares the term above and says that G^T·y and G^T· G belong to each other, (are not separable) and then takes furthermore the derivative, one gets the following result: 2 ·· (G^T· G). When one sets this term to zero since one is interested in the minimum and solves it after a one gets the result from Equation 30.

6.2.2 Normal equations

(30)

6.3 Singular-value decomposition (SVD)

A singular-value decomposition is one of the most widely used matrix operations in applied linear algebra.

6.3.1 Idea

Every Matrix G with the dimensions (N+ 1) × (m+ 1) can be decomposed as the triple product UDV ^T whereas U is an orthogonal (N +1) × (N +1) -matrix, D is a (N +1) × (m +1) diagonal matrix and V again is orthogonal with dimensions (m +1) × (m +1). When a matrix is orthogonal, the following applies: Q^TQ=QQ^T =I and Q^T =Q⁻¹. Due to that nice property, Equation 27 can be calculated according to Equation 32.

(31)

Ga_|=y_| becomes UDV^Ta_|=y_|. Which results in Equation 32

(32)

6.3.2 Uniform arguments and orthogonal polynomials

With uniform arguments x_i−x_j= (j−i)h for all i,j→=_t=0…N and orthogonal polynomials GG^T can be diagonalized and therefore the equation can be easier solved. An example can be found in subsubsection 6.3.6

(33)

with t= = =nCr(k,i)

p_k,N can now with g_k be put in the design matrix and the product G^TG will be a (m+ 1) × (m+ 1) Diagonal matrix. Afterwards a can be calculated with the knows formula G^TGa=G^Ty.

6.3.3 Calculation of the first terms for orthogonal polynomials:

From Equation 33 one knows that N⁽ⁱ⁾ and t⁽ⁱ⁾ are the following:

(34)

(35)

6.3.4 Exercise one, least square parabola

Compute a linear least-squares approximating parabola for the "second window" of five consecutive points (starting with x = 2) in the data.

Lets also say the following:

• m = deg = 2 (polys) {1;x;x²}
• N=4 (number of arguments actually 5 points)

Note: Normally: m≪N (The degree is smaller than the number of points)

Write down the design matrix and the system of normal equations.

(36)

With Equation 30 we get the system of normal equations as one can see in Equation 37:

(37)

Solve the linear system

When one wants to increase the stability of the matrix one can make a statistical normalization.

Compute the output y and the derivative (!) of the approximation at the central coordinate (x = 4)

6.3.5 Excercise three, Savitzky Golay filter

Apply the filter formulas developed in the exercise before for the data to compute approximately y for x= 3,…,8

6.3.6 Exercise four, orthogonal polynomials

Solve subsubsection 6.3.4 again by using the orthogonal polynomials _k=0,…,2
From Equation 33 and subsubsection 6.3.3 one knows the three basis functions:

Now we fist have to find a transformation. The transformation can be written in the follwoing way:

(38)

With Equation 30 we get the system of normal equations as one can see in Equation 39:

(39)

When we now calculate (G^TG)⁻¹ ·y one gets The result is therefore:

Which is the same as inparagraph 6.3.4.

6.3.7 Exercise five, singular value decomposition

Examine and compute a least-squares approximative quadratic parabola for the data

with respect to the basis functions in the following sense:

Compute the design matrix G and the normal matrix. Hint: The normal matrix here is diagonal!

(40)

(41)

Solve the system of normal equations and write down a formula for the approximating parabola.

divide first row by factor 5.

Also divide other rows by its factor.

What are the dimensions of the unitary matrices U,V , as well as the diagonal matrix D, in the singular value decomposition G=UDV^tr

• U=5x5
• D=5x3
• V=3x3

What are the entries (singular values) in the matrix D from above The singular values are the square-root of the non zero eigenvalues of G^T· G and therefore {;;}

Give three orthogonal basis polynomials (with respect to the data given) as formulas in the variable x. is orthogonal because G^T· G is diagonal.

6.4 Chebyshev polynomials

6.4.1 Idea

Approximate a continuous polynomial by a Chebyshev polynomial.

6.4.2 Definition

Chebyshev Polynomials are defined as T_n(x) = cos(n arccos(x)) with (n= 0,1,…) and (−1 ≤x≤ 1). Due to that, most data points are at the edge. The first polynomials can be found in Equation 42.

(42)

Further polynomials can be calculated with the recursion formula T_n+1(x) = 2xT_n(x) −T_n−1(n≥ 2) with the initial conditions T₁(x) =x,T₀(x) = 1.

6.4.3 Properties

• The maximal amplitude of a Chebyshev-Polynomials is and for a normalized one T_n+1(x)
• Amplitude: T_n(x) ∈ [−1,+1]
• Zero points: T_n(x) = 0 ⇔x= cos i= 0,1,…,n− 1 (also called chebyshev knots)
• T_n(x) =±1 ⇔x= cos(i= 0,1,…,n)

6.4.4 Usage

For Chebychev we use not g but T_x as function.

The matrix G^TG can be callculated accroding to Equation 43.

(43)

And results then in the follwing matrix

6.5 Continuous Chebyshev approximation

A weight function w(x) is often called the function inside the integral, in the case of the formula below w(x) =.

(44)

Equation 44 is true because the polynomial is orthogonal!

Example:Chebyshev continuous least-squares parabola on the interval [0,1] for y(t) =t³ First the interval [0,1] is transformed to [−1,1] by x= 2t− 1 This is shown with ??

In our case our new x is called t. Therefore x= 2t− 1 ⇒t= ⇒y(t) = Now we can use Equation 45 and get the following results:

(45)

6.6 Continuous Least-Square Legendre approximation

Legendre Polynomials: The Legendre polynomials are defined by the Rodriguez formula which can be seen in Equation 46

(46)

Where the first polynomials can be seen in Equation 47

(47)

Continuous Legendre least-squares approximation
If y(x) is function on [−1,1] which is absolutely square-integrable with respect to the weight function w(x) = 1 (x∈ [−1,1]) in the sense that ∫ ₋₁¹|y(x)|²dx<∞, then the continuous square-sum of residuals in Equation 48

(48)

Is minimal when the coefficients a have the value given in Equation 49. This equation also describes the resulting polynomial.

(49)

where m is the degree.

6.6.1 Legendre continuous least square parabola

Legendre continuous least-squares parabola on the interval t∈ [0,1] for y(t) =t³.

First one has to do a transformation to the interval x∈ [−1,1]. This results in the following: x= 2t− 1 ⇒t=, therefore y(x) =. With Equation 49 one gets then the following:

6.7 Mulit-variate least-square

Note: the degree of a multi-variate polynomial can be identified by adding up the degrees of the variables in each of the terms. It does not matter that there are different variables. The largest number is the degree.

• Input bi-variate (x,y) 2D : x
• Output scalar Z: 1D

• Degree 3.

  

• Degree 4.

  

Residuals: r_i=z_i−f (x_i) ⇒∑r _i² →min!

Basis Functions d=dim=2

Product

Statistical norm

6.7.1 Example one

Normally, a set of four 3d-points (x,y,z) is not contained in one single plane. But generally there is a plane coming close to the 3d-points in the sense of least-squares approximation.
The (x,y,z)-data in this problem is: A = (1,0,0), B = (0,1,0), C = (0,2,-1), D = (1,3,1)

(a)

Give a reasonable set of basis functions. Hint: A plane has total degree 1 (complete basis)

As a basis function, one can use:

(b)

Write down the design matrix according to a) and the normal equations

(c)

Solve the system of normal equations and give a functional formula for the approximating plane.

According to Equation 30 to following is true: _{normal equations}

6.7.2 Example three

Express x³,x⁴ as linear combinations of the Chebyshev polynomials T₀(x),T₁(x),T₂(x),T₃(x),T₄(x).

From Equation 42 one knows that T₃(x) = 4x³ − 3x⇒T₃(x) =x³ −x⇒T₃(x) +T₁(x) =x³


From Equation 42 one knows that T₄(x) = 8x⁴−8x²+1 ⇒T₄(x) =x⁴−x²+ ⇒T₄(x)+T₂(x) =x⁴− ⇒T₄(x) +T₂(x) +T₀ =x⁴

6.7.3 Example six

Express x³ as linear combinations of the Legendre polynomials P₀(x),P₁(x),P₂(x),P₃(x).
From Equation 47 one knows that P₃(x) =⇒P₃(x) =x³ −x⇒P₃(x) +P₁(x) =x³

6.7.4 Example seven

Compute continuously approximating least-squares lines for the model function y(t) =t² (0 ≥t≥ 1) by

(a)

Chebyshev approximation with the weight function w(x) = 1/ (−1 <x< 1)

Firstly one has to bring it into the correct range. One can do that with ??
Map the interval [a,b] onto the interval [c,d]

(50)

(b)

Legendre approximation with the weight function w(x) = 1

(c)

Estimate the maximum approximating errors in ab) by the coefficients of the Chebyshev polynomial T₂(x) and Legendre polynomial P₂(x) , respectively

7 Differentials, Taylor formulas and Jacobian

7.1 Differential

7.1.1 Definition

The purpose of differential is to measure error propagation. (How much is y(dependent variable) wrong when x(independent variable) is wrong by a certain amount). The differential df is the linear amount of change between a variable and a function, as it can be seen in Equation 51. Whereby ∆f is the hole amount of change, not only the linear one (The difference between two points). For small dx on can say ∆f ≈df and ≈.

(51)

7.2 Taylor

As one has seen before, ∆f ≈df for small dx and one has used only the linear part. To further improve the approximation one could not only use the linear part (first derivative) but also the squared (second derivative) and so on. Therefore, a function at a certain point can be approximated by it’s derivates at this point, which is called Taylor series approximation, which can be seen in Equation 52 for the one dimensional case and in Equation 53 for the multidimensional case.

(52)

The vector field of of the partial derivative is called gradient of f and is denoted with grad(f ) or ∇f (x). Therefore ∆f ≈∇f (x) ·h. Equation 53 shows the taylor series approximation for the multi-indices α= and |α|=α₁ +…+α_n.

(53)

The remainder terms R_α are absolutely bounded by max_x∈S with |α|=N+ 1 and S=x₀ + is a n-dimensional ’rectangle’ with center x₀. (This idea can be used afterwards by the Jacobian matrix and the determinant)
The formulas for up to order four can be found in Equation 54

(54)

7.2.1 Example

The bivariate symmetric (!) function f (x,y) =e⁻ has to be approximated by a bivariate Taylor polynomial of order 4 around (0,0) by

1.

evaluating and using Equation 53 for the partial derivatives.

2.

multiplying the univariate Taylor series for the exponential function,

3.

substituting into the univariate Taylor series for the exponential function

4.

What value do you expect for the error limit: lim_∥h∥→0

7.3 Jacobian matrix and determinant

In the Taylor series approximation for multi-indices, one has seen that ∆f ≈∇f (x) ·h for small dx (see also Equation 53). When one writes all those gradients in one matrix, one gets the so-called Jacobian matrix. Which actually tells us the same as mentioned in subsubsection 7.1.1, but just for a multivariable problem. It tells us how y₁,y₂,… (dependent variable) changes when, x₁,x₂,… (independent variable) changes. Furthermore it has the nice property that the determinant of this matrix also describes how the volume changes when one changes a certain variable.Therefore the Matrix can be used to analyse the error propagation or make volume calculation in a different coordinate system.
Below one can find some definitions.

(55)

(56)

(57)

In the special case that n = m the Jacobian matrix is a square matrix and thus has a determinant (called Jacobian determinant):

(58)

(59)

7.3.1 Estimating navigation error by inversion of Jacobian determinant

Lets assume one measures an angle (α,γ) and want the position (x,y). Let’s call the transformation from (α,γ) to (x,y) T and from (x,y) to (α,γ) T⁻¹

To solve this problem one can follow the following approach:

1.

Description of the problem

2.

Transformation (equalities) from cosine-theorems (generalized Pythagoreans)

3.

Implicit differentiation

4.

Computation of the Jacobian matrix J_T⁻¹(x,y)

5.

Computation of the Jacobian determinant

6.

Elimination of angles

7.

Computing expressed in position coordinates x,y

7.3.2 Example three

The formulas x=r cos(φ), y=r sin(φ) (0 <r, 0 ≤φ< 2π) define the coordinate transform T from polar coordinates (r,ϕ) to Cartesian coordinates (x,y) in the punctured plane R²' {(0,0)}.
Compute the Jacobian matrix J_T(r,φ) and its Jacobian (determinant) detJ_T(r,φ) as well as the Jacobians (determinants) detJ_T(x,y),detJ_T⁻¹(r,φ) and detJ_T⁻¹(x,y)

First of all one has to write down the transformations

Once one has done that one can calculate the Jacobian Matrix:

and afterwards the determinant

when one wants to express r in x and y one knows that r= therefore

and

7.3.3 Example one

Elliptical coordinates (σ,τ) (σ> 1,−1 <τ< 1) for a = 1 are connected to Cartesian coordinates (x, y) through the transforming formulas:

Compute the Jacobian matrix

Express the Jacobian determinant

Express the Jacobian determinant in x, y, since the denumerator is exactly y, we know it already. Furthermore the following is true: =σ² −τ²

8 Ordinary differential equations

8.1 Definition

An ODE is a differential equation where it’s derivatives belong to only one variable. Furthermore it is possible that an ODE can not be solved explicitly, due to that in this chapter it will be investigated how they can be solved numerically.
A first order differential value problem can be expressed like the one in Equation 60

(60)

8.2 Explicit methods

One way to solve the problem numerically is by tailor series approximation, as it can be seen in Equation 61. (Remember that the tailor series is evaluated from one single point). The idea is that one creates a tailor series approximation up to order p at the starting point and moves then with a step size h to the new position one got from the tailor series approximation and does there the same until one reaches the destination.

(61)

The calculations up to order three (p=3) are given in Equation 62:

(62)

8.2.1 Euler method

The Euler method is a special case of the explicit methods where p= 1 and h=const. The Formulas to calculate it can be found in Equation 63.

(63)

8.2.2 Error Calculation

• Global Error:

  (64)

• Local Error:

  (65)

• Local slope Error:

  (66)

8.2.3 Example

Solve the initial value problem y^′=xy^1/3 with y(1) = 1 numerically by the method of Taylor with order p = 4 and fixed step-size h = 0.1 for the x-values 1.1 and 1.2 (two steps). All final (!) results should be rounded to the 10th digit. Furthermore, compute the local error (slope) as well as the global error for the two steps. Note that the exact solution of the equation is given by Equation 67.

(67)

To find out x at 1.1 one firstly needs to calculate the it’s tailor series approximation at the starting point.

With Equation 61 one can then write down the following:

Then one does the same at the new location on got from the previous result.

The global error is defined by Equation 64 and therefore has the following result:

The local error is defined by Equation 66 and therefore has the following result:

and for n=1 one has the following result:

8.3 Explicit Runge-Kutta Methods

Since the calculation of the tailor approximation series is quite tedious. One came up with another method, which is recursive and results in the same as the tailor series approximation (big advantage is that no computations of derivatives up to order p are needed). In Equation 68 one can find the Common Runge-Kutta method of order 4. As one can see there one has seven variables, when one now sets this equal with the tailor series approximation all those variables are defined in the end as it can be seen in Equation 69.

(68)

(69)

For the classic runge-kutta method one can also write Equation 70 which results in Equation 71.

(70)

(71)

The errors are calculated nearly the same way as in the previous section. Equation 74.

• Global Error:

  (72)

• Local Error:

  (73)

  When there is only one step calculated, the local error must be the same as the global error.

• Local slope Error:

  (74)

8.3.1 Example

Solve the initial value problem y^′=xy^1/3 y(1) = 1 numerically by the classical Runge-Kutta method of order p = 4 and fixed step-size h = 0.1 for the x-values 1.1 und 1.2 (two steps). All final (!) results should be rounded to the 10th digit. The exact solution of the equation is y=^3/2.
From Equation 71 one knows that:

As one can see, the result is nearly the same as in the previous Example. The local slope error can be calculated according to Equation 74.

The local error can be calculated according to Equation 73 and results in the following:

8.4 Butcher tableau

The Butcher tableau is mainly a mnemonic device (Gedächtnissstütze) to remember the coefficients. The tableau must fullfill the condition listed in Equation 75

(75)

Where the variables mean the folwoing:

• s= number of stages
• c= ’inputs’ Subdivides inputs
• b= weights
• a = subdiv ’y’ ’output’

The general tableau can be seen in Equation 76.

(76)

A step from x_n to x_n+1 =x_n+h_n(n= 0,1,… ) in the general Runge-Kutta method is defined by Equation 77 where the values a,b,c can be read from the butcher tableau.

(77)

Some examples for butcher tableaus can be found in Table 1

The simplest adaptive Runge–Kutta method involves combining Heun’s method, which is order 2, with the Euler method, which is order 1 (also called Heun-Euler 2(1)). Its extended Butcher Tableau can be seen in Equation 78.

(78)

8.5 Step-size adaption

8.5.1 Idea

The Idea is to automatically adapt the step-size h. Due to that one needs a new way to define the approximation of the error, which can be done with an Accuracy Goal (ac), which defines how many decimal places (Nachkommastellen) are correct and a precision Goal (pg), which represents the significant digits of the result. The two parameters are considered in the tolerance parameter, ε which can be found in Equation 79.

(79)

Furthermore, one needs a second approximation for the error calculation, with the order ˆ p. The first order approximation p is needed to calculate the step size. Mostly ˆ p=p− 1. Due to that the butcher tableau is extended by a row (b values) as it can be seen in Table 2.

0	0	0		0	0
c 2	a2,1	0		0	0
				0
cs−1	as−1,1	as−1,2		0	0
cs	as,1	as,2		as,s−1	0
	b1	b2		bs−1	bs
	ˆ b1	ˆ b2		ˆ bs−1	ˆ bs

The local error can be calculated according to Equation 80.

(80)

Below one can find again a short description of the variables:

• ε_a=dy: absolute error ε_a= 10^−ag
  Ex. ag= 4 ⇒ 0.001289
• ε_r=: relative error ε_r= 10^−pg
  Ex. pg= 4 ⇒ 0.00_{signific. digit}

To know if the step size is good or not one calculates . When the current step > 1, then the estimation of h_n was too optimistic and the step must be repeated with a smaller step size. One also says the current step is Rejected. Otherwise when ≤ 1 the step size is ok and one can proceed. Updating the step size is done according to Equation 81.

(81)

ε=ε_a+ε_r with ˜ p= min(p, ˆp) + 1 (order of the primary method)

8.5.2 Stability of explicit methods

The global relative error must not diverge, which means it must be limited. Since it is difficult to make a statement about the analysed ODE one uses benchmark equations. One of the most commonly used ones is the Dahlquist model, which can be seen in Equation 82

(82)

The solution of this equation is y=e^ℜ{A}x, which is an oscillation with a exponential amplitude e^ℜ{A}x and frequency ℑ{A}. The

Example Euler

The exact solution is: Y(x) =e^−λx Consider Euler’s method: Solation will go to zero iff

Stability for heun-method (rk-2)
From Equation 77 one knows that y_k+1 =y_k+h· ( ·k₁ + ·k₂)

(83)

From that one can somehow derive three cases which are listed below

3 Cases

• ℜ{A}< 0: The Amplitude of the ODE gets smaller exponentially. The stability conditions says that the approximated values y_k(k = 0,1,…) must be exponentially damped too. Therefore and from y_k+1 ∝|F(z)| the following must hold |F(z)|< 1.
• ℜ{A}= 0 : The Amplitude of the ODE is constant. The stability conditions says that the approximated values y_k(k = 0,1,…) are constant too. Therefore and from y_k+1 ∝|F(z)| the following must hold |F(z)|= 1.
• ℜ{A}> 0 : The Amplitude of the ODE gets larger exponentially. The stability conditions says that the approximated values y_k(k = 0,1,…) must be exponentially growing too. Therefore and from y_k+1 ∝|F(z)| the following must hold |F(z)|> 1.

The stability condition for case one can know be calculated as the following (ℜ{A}< 0) and 1 >|F(z)|=⇒−2 <hℜ{A}< 0. since x_1,2 == =−1 ± 1 and therefore the values must be between 0 and -2. The stability polynomial in this exercise was F(z) = 1 +z+ (z=hA).

Recursive Formulas
For the stability polynomials in the form: F(z) = 1 +b₁k₁(z) +…+b_sk_s(z) recursive equations exist as can be seen in Equation 84.

(84)

8.5.3 Exercise adaptive step size

Solve the initial value problem φ^′=c(1 −εcosφ)² φ(0) = 0 with c= 1 and ϵ= 0.25 numerically by applying the Heun-Euler 2(1) embedded adaptive method with classical step-size control until 3 proceeding steps are executed. The initial step-size equals 0.001, the accuracy goal (ag) 4 and the precision goal is 4 , either.
Create a table listing values for , state ) containing at least three preceding steps.
According to the exercise, we know the following:

1.

The position at the beginning: x= 0, y= 0

2.

The step size: h= 0.001

3.

The local error can be calculated according to Equation 80. Which says that e_k=h_n∑ _j=1^sk_j. From Equation 78 one knows b₁ =, ˆ b₁ = 1, b₂ = and ˆ b₂ = 0. Furthermore, one knows from Equation 77 and Equation 78 that k₁ =f= 1 ·² = and k₂ =f= 1 ·² =0.5625. Therefore e_k =h_n· (b₁ − ˆ b₁) ·k₁ + (b₂ − ˆ b₂) ·k₂ =_hn· (_{b1− ˆb1} · +_{(b2− ˆb2)} ·0.5625) =2.966 · 10⁻¹¹

4.

Since p= 2 (order) and ε can be calculated according to Equation 79 which says: ε=ε_a+|y|ε_r = 10^−ag+|y|10^−pg= 10⁻⁴ + 0 · 10⁻⁴ = 10⁻⁴. Therefore ⁻ =⁻ = 1.836 · 10³

5.

With Equation 81 one can the finally calculate the new step size which is: h_n+1 = h_n⁻ = 0.001 ·⁻ = 1.836

6.

the new y value can be calculated according to Equation 77 which means for the current scheme 0.001 · ( ·k₁ + ·k₂) = 0.001 · ( · + ·0.5625) = 0.0005625

x	y	hn	ek	−	hn+1	state
0	0	0.001	2.966 · 10−11	1.836 · 103	1.836	Proceed
0.001	0.0005625	1.836	0.18169	0.02347	0.043	Reject

8.5.4 Exercise Stability polynomial

Using Theorem 1.3 in the script (p. 29) gain the A-stability polynomial F₁(z) = 1 +z+ + for the embedded adaptive method SS3(2) with Butcher tableau.

From Equation 84 one knows that a stability polynomial is of the form F(z) = 1+b₁k₁(z)+b₂k₂(z)+b₃k₃(z)+0k₄(z), whereas:

8.5.5 Stiffness

The stiffnes is dependent on:

• ODE
• step size h
• Method (for example RK-4)

Stiffness Detection
Heun euler does not work when we have stiffness situations. Condition: c_s−1 =c_s= 1. Through testing of to the absolute borders of the stability region the stiffness can be detected. Stiffness is present when |h˜λ| is outside or at the border of the stability condition.

Example for explicit runge-kutta

˜λ is an estimation for f_y=f (x,y) for example y^′=x⁴ − 25y⁴ =f (x,y) and takes the role of ℜ{A} for the stability analysis.

8.5.6 Exercise stiffness detection test

Given the differential equation y^′=_{f (x,y)} y(0) = 4 carry out a stiffness detection test using the A-stability region and the partial derivative f_y at the initial values and step size h= 1. The method is defined by the Butcher tableau below.

The exercise can be solved in two steps:

1.

Get stability region (Equation 84):
From Equation 84 one knows that a stability polynomial is of the form F(z) = 1 +b₁k₁(z) +b₂k₂(z) + b₃k₃(z) + 0k₄(z), whereas:

To calculate the stability region one hast to set the stability polynomial to zero and therefore one gets: x_1,2 == =−1 ± 1 and therefore the boundary is -2 and 0.

2.

Check stiffness (see subsubsection 8.5.5)
By calculating the partial derivative of − after y one obtains the following −1 · · · 2 ·y inserting x= 0 and y= 4 results in −1 ⇒A=−1 and h·A=−1 which is inside the A-stability constraint −2 <Ah< 0 and therefore no stiffness is detected.

8.5.7 Van der Pol second-order differential equation

Solve the van-der-Pol ODE-system z(0) = 1,v(0) =−1, µ= 0.2 numerically by applying the Heun-Euler 2(1) embedded adaptive method with classical step-size control until 3 proceeding steps are executed. The initial step-size equals 0.001, the accuracy goal (ag) 1 and the precision goal is 2. Create a table listing values for state) containing at least three proceeding steps.

According to the exercise, we know the following:

1.

The position at the beginning: x= 0, z= 0,v=−1

2.

The step size: h= 0.001

3.

The local error can be calculated according to Equation 80. Which says that e_k=h_n∑ _j=1^sk_j. From Equation 78 one knows b₁ =, ˆ b₁ = 1, b₂ = and ˆ b₂ = 0. Furthermore, one knows from Equation 77 and Equation 78 that k₁ = f== = and

Therefore e_k= 0.001 · (− ·+ ·) =

4.

Sine p = 2 (order) and ε can be calculated according to Equation 79 which says: ε=ε_a+|y|ε_r= 10^−ag+|y|10^−pg= 10⁻¹ +· 10⁻² =. Therefore ⁻ = = 469.04157598235

5.

With Equation 81 one can the finally calculate the new step size which is: h_n+1 = h_n⁻ = 0.001 · 469.04157598235 =0.46904157598235

6.

the new y value can be calculated according to Equation 77 which means for the current scheme 0.001 · ( ·k₁ + ·k₂) = 0.001 · ( ·+ ·) =

x	y	hn	ek	−	hn+1	state
0		0.001		469.04157598235	0.46904	Proceed
0.001		0.46904157598235				Reject

(85)

(86)

(87)

(88)

(89)

Where

Common angles

Degrees	0◦	30◦	45◦	60◦	90◦
Radians	0
sin 𝜃	0				1
cos 𝜃	1				0
tan 𝜃	0		1

Reciprocal functions

cotx	=
cscx	=
secx	=

Even/odd

Pythagorean identities

Cofunction identities

sin( −x)	= cosx
cos( −x)	= sinx
tan( −x)	= cotx
cot( −x)	= tanx
sec( −x)	= cscx
csc( −x)	= secx

Sum and difference of angles

sin(x+y)	= sinx cosy+ cosx siny
sin(x−y)	= sinx cosy− cosx siny
cos(x+y)	= cosx cosy− sinx siny
cos(x−y)	= cosx cosy+ sinx siny
tan(x+y)	=
tan(x−y)	=

Double angles

sin(2x)	= 2sinx cosx
cos(2x)	= cos2x− sin2x
	= 2cos2x− 1
	= 1 − 2sin2x
tan(2x)	=

Half angles

sin	=±
cos	=±
tan	=
	=

Power reducing formulas

sin2x	=
cos2x	=
tan2x	=

Product to sum

sinx siny	=cos(x−y) − cos(x+y)
cosx cosy	=cos(x−y) + cos(x+y)
sinx cosy	=sin(x+y) + sin(x−y)
tanx tany	=
tanx coty	=

Sum to product

sinx+ siny	= 2sincos
sinx− siny	= 2cossin
cosx+ cosy	= 2coscos
cosx− cosy	=−2sinsin
tanx+ tany	=
tanx− tany	=

9 Fromulas

9.1 Differentation Formulas

9.2 Integration Formulas

9.3 Table of Indefinite Integrals

9.3.1 Basic Functions

9.3.2 Products of e^x and cox x and sin x

8.

∫ e^axsin(bx)dx=e^ax(a sin(bx) −b cos(bx)) +C

9.

∫ e^axcos(bx)dx=e^ax(a cos(bx) +b sin(bx)) +C

10.

∫ sin(ax)sin(bx)dx=(a cos(ax)sin(bx) −b sin(ax)cos(bx)) +C,a≠b

11.

∫ cos(ax)cos(bx)dx=(b cos(ax)sin(bx) −a sin(ax)cos(bx)) +C,a≠b

12.

∫ sin(ax)cos(bx)dx=(b sin(ax)sin(bx) +a cos(ax)cos(bx)) +C,a≠b

9.3.3 Product of Polynomial p(x) with lnx,e^x,cosx,sinx

13.

∫ xⁿlnx dx=xⁿ⁺¹ lnx−xⁿ⁺¹ +C,n≠− 1

14.

∫ p(x)e^ax dx=p(x)e^ax−∫ p^′(x)e^ax dx=p(x)e^ax−p^′(x)e^ax+p^′′(x)e^ax−…
(Sign alternate: +−+−+−… )

15.

∫ p(x)sin(ax)dx=−p(x)cos(ax) +∫ p^′(x)cos(ax)dx =−p(x)cos(ax) +p^′(x)sin(ax) +p^′′(x)cos(ax) −…
(Sign alternate in pairs after first term: −++−++… )

16.

∫ p(x)cos(ax)dx=p(x)sin(ax) −∫ p^′(x)sin(ax)dx =p(x)sin(ax) +p^′(x)cos(ax) −p^′′(x)sin(ax) −…
(Signs alternate in pairs: ++−−++−… )

9.4 Taylor Polynomial/Series

Development of f around a

in which k! = 1 · 2 · 3 ·…·k.

9.4.1 Important Taylor Series

The last three only converge for |x|< 1.

9.5 Determinant

9.5.1 Sarrus

9.6 Matrix

9.6.1 Transpose

(90)

9.6.2 Multiplication